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ABSTRACT 

The linearly-polarized solar limb spectrum that is produced by scattering pro- 
cesses contains a wealth of information on the physical conditions and magnetic 
fields of the solar outer atmosphere, but the modeling of many of its strongest 
spectral lines requires solving an involved non-LTE radiative transfer problem 
accounting for partial redistribution (PRD) effects. Fast radiative transfer meth- 
ods for the numerical solution of PRD problems are also needed for a proper 
treatment of hydrogen lines when aiming at realistic time- dependent magneto- 
hydrodynamic simulations of the solar chromosphere. Here we show how the 
two-level atom PRD problem with and without polarization can be solved accu- 
rately and efficiently via the application of highly convergent iterative schemes 
based on the Gauss-Seidel (GS) and Successive Overrelaxation (SOR) radiative 
transfer methods that had been previously developed for the complete redistri- 
bution (CRD) case. Of particular interest is the Symmetric SOR method, which 
allows us to reach the fully converged solution with an order of magnitude of 
improvement in the total computational time with respect to the Jacobi-based 
local ALI (Accelerated Lambda Iteration) method. 

Subject headings: line : profiles - methods: numerical - polarization - radiative 
transfer - scattering - stars : atmospheres - Sun: atmosphere 



1. Introduction 



The phenomenon of scattering in a spectral line is a complicated physical proce ss where 
partial correlations between the incoming and outgoing photons can occur (e.g., iMihalas 
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19781 ; ICannonl Il985t lOxenius &: Simonneaul Il994l ). This happens when the shape of the 
incident spectrum that populates the upper level via radiative absorptions is not flat over 
the line. These partial redistribution (PRD) effects tend to be more important in strong 
lines, such as the solar Mg II and Ca II resonance lines and Lyman a. In particular, the 
wings of the intensity profiles of these lines are strongly affected by PRD effects, especially 
concerning observations close to the edge of the solar disk. 

Only a small number of solar spectral lines show conspicuous PRD signatures in their 
emergent intensity profiles. However, a substantially larger fraction show clear hints of PRD 
effects in the fractional linear polarization Q/I profiles that result from scattering processes 
in quiet regions of the solar atmosphere (e.g., see the classification proposed by Belluzzi & 
Landi Degl'Innocenti 2009 of the various Q/I shapes found in the linearly-polarized solar 
limb spectrum observed by Stenflo & Keller 1997 and by Gandorfer 2000, 2002, 2005). To 
achieve a rigorous modeling of the weak polarization signals that constitute this so-called 
Second Solar Spectrum is very important, mainly because it contain s valuable informatio n on 
the magnetism of the extended solar atmosphere (e.g., the review by lTrujillo Buenoll2009l ). To 
this end it is crucial to solve accurately and efficiently the non-LTE (Local Thermodynamic 
Equilibrium) radiative transfer problem of resonance line polarization taking into account 
PRD effects. The present paper represents a contribution towards this goal. 

Fas t iterat ive methods based on operator splitting were introduced to astrophysics by 



Cannon! (119731 ) for unpolarized radiative transfer with complete redi stribution (CRD) in 
scatte ring. Extensio ns of this type of met hods to PRD were done by IVardavas Cannon 
(119761 ) . and later by IScharmerl ( 119831 ) and lUitenbroek! (120011 ) . These methods ar e widely 
know n today as Accelerated Lambda Iteration (ALI) methods (e.g., the review by iHubeny 
20031 ). An optimum choice for the approximate lambda operator is the diagonal of the 
true la mbda operator, which was introduced in the seminal paper by lOlson. Auer fc Buchler 
( 119861 ). This Jacobi based me t hod f or solving the two- level atom problem with CRD was 
generalized by lPaletou fc Auerl ( 119951 . hereafter PA95) to unpolarized PRD radiative transfer. 



Superior radiative transfer methods ba sed on Gauss-Seidel (GS) and succe s sive o verre- 
laxation (SOR) iteration were developed by lTrujillo Bueno &: Fabiani Bendichol (119951 here- 
after TF95). The convergence rate of these iterative schemes is equivalent to that corre- 
sponding to upper or lower triangular approximate lambda operators, but without the need 
of constructing and inverting such operators. Therefore, the computing time per iteration 
is similar to that of the Jacobi scheme or local ALI method, but the number of itera- 
tions needed to reach convergence is an order of magnitude smaller. In their paper, TF95 
suggested the strategy to generalize their GS-b ased methods to the multilevel atom case. 
Fabiani Bendicho. Trujillo Bueno fc Auerl ( 119971 ) implemented such a Multilevel GAuss Sei- 
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del method (MUG A; see also Fabiani Bendicho & Trujillo Bueno 1999 and Asensio Ramos Sz 
Trujillo Bueno 2006 for its generalization to 3D and spherical geometries) and combined it 
with a non-linear multigrid iterative scheme to produce multilevel radiative transfer programs 
whose convergence rates are insensitive to the grid size. Here we present the generalization 
of the GS and SOR radiative transfer methods of TF95 to the two-level atom PRD problem, 
with and without scattering polarization. 

An alternative iterative scheme for solving radia tive transfer problems is the precon- 
ditioned bi-co njugate gradient method (ICastorl I2004J) , which ha s been recently applied to 
plane-parallel ( Paletou fc Anterrieu 20091 Nagendra et al. 2009 ) and spherical geometries 
( Anusha et al. 20091 ). Its rate of conv ergence is similar to that of the optimal Symmetric 



SOR method of TF95 (see ICastorll2004l ) . but an efficient generalization of the preconditioned 



bi-conjugate gradient method to the multilevel atom case is presently an unsolved problem. 

A suitab le generalization of the local ALT method to t he Zeeman line transfer problem 
was done by iTrujillo Bueno fc Landi Degl'Innocentil ( 119961 ) . Extensions of the Jacobi, GS , 

and SOR sche mes to scattering polarization were carried out by Faurobert-Scholl et al.l ( 1997 . 

CRD Jacobi) JPaletou fc Faurobert-Scholll ( 119971 . PRD Jacobi), and ITrujillo Bueno fc Manso Sainz 
( 119991 . CRD Jacobi, GS and SOR). All these iterative schemes for non-LTE radiative transfer 
were generalized to solve CRD multi-level scattering polarization problems in the presence of 
a magnetic field, including the possibility of atomic polarization in all levels (Trujillo Bueno 
1999; Manso Sainz & Trujillo Bueno 2003, see also Trujillo Bueno 2003). However only the 
Jacobi iterative scheme has been applied to solve the two-level atom polarized PRD radiative 
transfer problem in the presence of an external magnetic field (Nagendra et al. 1999; Fluri 
et al. 2003; Sampoorna et al. 2008, see also the reviews by Nagendra 2003; Nagendra & 
Sampoorna 2009). Here we extend the GS and SOR iterative methods to solve polarized 
PRD problems in the absence or in the presence of magnetic fields which do not break the 
axial symmetry of the problem (e.g., the micro-turbulent field case). 

The accuracy of any iterative method for a given depth grid resolution is determined b y 



the truncation error or the true error (see lAuer. Fabiani Bendicho fc Trujillo Bueno! 119941 ) . 



So far the study of the true error is limited to only CRD problems (e.g., Auer et al. 1994; 
TF95; Faurobert-Scholl et al. 1997; Trujillo Bueno & Manso Sainz 1999; Chevallier et al. 
2003). Hence, in this paper we discuss certain aspects of the true error for PRD problems. 

For our stu dy we consider all the three angle-averaged (AA) redistribution functions of 



Hummer! (119621 ). namely #1,11,111, aa, an d the linear combination of -Ri^aa and Ru^aa- We 
recall that physically (1) -R^aa represents the case of infinitely sharp lower and upper levels 
(or pure Doppler redistribution in the laboratory frame); (2) Rx^aa represents the case of 
infinitely sharp lower level and radiatively broadened upper level (coherent scattering in the 
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atomic frame); (3) -Rhi,aa represents the case of infinitely sharp lower level and radiatively 
as well as collisionally broadened upper level (CRD in the atomic frame). 

The logical structure of this paper is the following: §§ 2, 3, and 4 are devoted to 
unpolarized PRD radiative transfer. We first recall the basic equations, the Jacobi scheme 
used by PA95 for i?ii,AA redistribution, and then present briefly the extension of this scheme 
to -Ri,iii,aa redistributions. Next, we present the generalization of the GS and SOR schemes 
of TF95 to PRD. A detailed study of the true error for all the three iterative schemes with 
AA redistribution functions is conducted in § 4. §§ 5, 6, and 7 are devoted to polarized 
PRD radiative transfer. In § 5 we present the basic equations of polarized radiative transfer. 
Our generalization of the Jacobi, GS, and SOR schemes of § 3 to scattering polarization is 
presented in § 6. A detailed study of the true error for polarized PRD case is given in § 7. 
Concluding remarks are given in § 8. 



2. Unpolarized PRD radiative transfer 

We consider the case of scattering on a two-level atom with a backgr ound cont i nuum 



The scattering mechanism is described by the AA redistribution functions of lHummerl (119621 ). 
Furthermore, we approximate the stellar atmosphere by a one- dimensional plane parallel slab 
of total optical thickness T. Under these assumptions, the scalar radiative transfer equation 
is given by 

j^-M 7 ") = J ^( r ) - 3 x (t), (1) 

where I x ^{t) is the specific intensity, x is the frequency from line center in units of the Doppler 
width, n = cos#, with 9 being the angle between the ray and the atmospheric normal, and 
the optical depth r is defined by dr = — {xi4>x + Xc)dz//i, with <p x the normalized Voigt 
profile function, \ c and xi t ne continuum and line opacities, and z the distance along the 
normal to the atmosphere. Hereafter, we omit the r dependence of the intensity and source 
function for notational simplicity. The monochromatic source function is given by 

_ (f> x Si x + rB 

where r = Xc/xu an d B is the Planck function at the line frequency. The line source function 
is given by 

S lx = (l-e)J x + eB, (3) 

where e is the collisional destruction probability. The PRD scattering integral or mean PRD 
intensity is given by 

Jx= gL'Jx'dx', (4) 
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with the mean intensity 



J r 



(5) 



In Equation 



9xx' 



Rk.\\(x,x' ) /(I) T , w here RkAA with k — I, II, and III are the AA 



Hummer 
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MihalasI f 


1978 


) and 


Heinzel 


(1981) 



3. Iterative methods for unpolarized PRD radiative transfer 



We write the formal solution of the radiative transfer equation ([I]) as 



Xfl 1 



(6) 



where T Xfl gives the transmitted specific intensity due to the incident radiation at the bound- 
ary and A x/1 is a iV x N operator whose elements depend on the optical distances between 
the grid points, with N being the number of spatial grid points. 

A suitable formal soluti on method for the numerical solution of Equation (ED) is th e 
short-characteristics method (IKunasz fc Auerlll988l ; lAuer fc Paletoulll994l ; lAuer et al.lll994l ). 
This method is based on the assumption that the variation of the source function with the 
optical depth along the ray under consideration is a parabola between three consecutive grid 
points. Thus, if M represents an upwind point, O represents the point of interest and P the 
downwind point, then the intensity at point O is given by 



lxn,0 — Ixti,M e + ^x,m(^)S X: M 

+ ^x,o(l X ) S x,0 + ^x,p{P>)S x ,p, 



(7) 



where At m is the optical distance on segment MO, ^ x ,m,o,p are functions of the optical 
distances between O and M and between O and P, and S Xt M,o,p are source function values 
at the points M, O, and P, respectively. However, at the boundaries a linear interpolation 
for the source function along the points M and O is used for the rays going out of such 
boundaries. 



Following iTrujillo Buend (120031 ) we write the mean PRD intensity at the spatial grid 
point V as 



J dxx'i^x' ,ilS x > ; i + ■ • • + ^-x' ,ii-l^x' 

+A x ijNS xltN ]dx' + T X}i , 
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where T X j is given by Equation (j3J) but with J x > replaced by T x > (which is given by Equa- 
tion ([5]) with I x i^ replaced by T x i^). For each frequency x, A Xiii > is obtained by integrating 
A^ over all the directions \i of the incoming and outgoing radiation beams of the angular 
quadrature chosen for the numerical integration. Furthermore, a, b and c are simply symbols 
that will be useful to indicate weather we choose the 'old' or the 'new' values of the source 
function. In the following subsections we, first, briefly recall the Jacobi iterative method, 
and then present the GS and SOR iterative schemes. 



3.1. Jacobi iterative scheme 

Let us recall first the Jacobi iterative scheme presented in PA95. This scheme is obtained 
by choosing in Equation (JSJ) a = c = old, but b = new, which gives : 



Jx,i — J x ]t + J g xx '^x',iiPx' SSi x i ti dx' , 



(9) 



where we have used Equation In the above equation p x = 4> x /{4> x + r) and SSi Xji = 
SfxJ ~ ^ix\- U sm g Equation (jUJ) in Equation we obtain the following expressions for the 
line source function corrections : 

-(l -<)/<£, W^d.- 

= (1 " e)J% + & - SgJ . (10) 
Following PA95 we define a frequency dependent residual as 

[l-ep^ + eB-S^. (11) 



X,l 



Thus at each depth point we have to solve N x linear equations, where N x is the number 
of frequency points. The simplest (but numerically expensive) way to find the solution of 
the system of linear equations (flOl) is by matrix inversion as follows : 

6S = A~ 1 r, (12) 

where at each depth point i, r is a vector of length N x , and the matrix A is of dimension 
N x x N x , and its elements are given by 

A mn = 5 mn — (1 — e)gm n A n ,uPn i m i n = 1; ■ ■ ■ > A x . (13) 

Here 5 mn is the Kronecker's symbol, g^ n are the redistribution weights, and the indices m 
and n refer respectively to the discretized values of x and x' . Note that for isothermal slabs 
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the matrix A can be computed only once and can be inverted and stored. This method was 
referred to as the Frequency-by- Frequency (FBF) method by PA95, which was developed by 
the authors for type II redistribution. It is easy to note that the same method can be easily 
applied to the type I and type III redistributions and a linear combination of type II and 
type III redistributions without any further effort. 

The above FBF method involves the inversion of a matrix, which can be huge for realistic 
problems. For this reason, PA95 proposed a faster but equally robust method for the case 
of a type II redistribution function. In the following sub-section we briefly discuss this more 
efficient method, which is presented in more detail in PA95. 



3.1.1. CRD-CS or Core-Wing method for type II redistribution 
It is well kn own that q^,, b ehaves like CRD in the line core and like coherent scattering 



in the wings (see lMihalaall978l ). Taking advantage of this fact, to reduce the computational 
cost involved in the calculation of 5Si x j, one can introduce a core- wing approximation to the 
true redistribution function g^ x ,, namely 

ii _ / <f>*, for x,x' < x c , 

5(x — x'), for x' > x c . 



1XX 



Here x c is called the separation frequency that distinguishes between the line core and the 
wing. PA95 showed that x c = 3.5 Doppler widths is a physically reasonable choice and, 
hence, we adopt the same in this paper. 

Substituting Equation (Tfjj) in Equation (ITUl) . the equation for 5Si Xt i takes the simpler 

form 

6S lxi = + p ~ ^ , (15 ) 
l-a x (l- e)p x A Xtii 

where a x are the splitting coefficients that allow for a smooth transition between the core 
and the wing. In the core a x = and in the wing a x = g\\. The frequency independent core 
integral ATj is given by 

ATj = (l-e) / 4> x 'p x 'K xlyi i 8Si X ! yi dx'. (16) 



To evaluate ATj, we consider only the core frequencies (i.e., a x = 0) in Equation ( fl5l) and 
then apply the operator (1 — e) J core 4>x'Px'^-x',u[]dx' . The resulting scalar equation for ATj 
can be easily solved to obtain ATj (see PA95 for more details). From Equation ffT5]) we see 
that the wing frequencies drop out for x < x c , and in the wings the term multiplying (1 — a x ) 
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appears only as a frequency independent quantity. Thus, it is possible to find all the 5Si x< i 
values from a simple scalar equation, thereby completely avoiding the solution of a system 
of equations irrespective of the number of frequency grid points. 



3.1.2. Extending the CRD-CS method to type I and III redistribution 



The extensio n of the CRD-CS method to type III redistribution has been given in 
Fluri et al.l ( 120031 ). To this end, the type III redistribution function is approximated by 
assuming CRD in the core and by setting it to zero in the wings. This is justified because 
the t ype III redistribution function does not show coherent peaks in the wings (see iMihalas 
19781 ). However, we find that for the pure type III redistribution case we have to approximate 
-Riii,aa by CRD throughout the frequency bandwidth to compute SSi Xt i. Setting it to zero in 
the wings leads to convergence problems. In the case of a linear combination of -Rh,aa and 
i?ni,AA (see § 4.4, Equation (1211 ). we find that as long as elastic collisions are small we can 
approximate Rj^aa by CRD-CS and -Rulaa by CRD in the line core and set it to zero in the 
wings. However, when elastic collisions are large we can approximate Rh,aa by CRD-CS, 
but -Riii,aa should be approximated by CRD throughout the line profile, otherwise we have 
convergence problems. 

For soling the type I redistribution problem we approximate Ri t AA by CRD in the core 
and set it to zero in the wings for the computation of the SSi Xt i corrections. This is justified 
as Rj t AA is a pure Doppl er redistribution function and does not show coherent peaks in the 



wings (see lMihalaslll978l ). 



3.2. Gauss-Seidel and SOR iterative schemes 

The radiative transfer methods based on GS and SOR iterations were developed by 
TF95 for the CRD case. In this section we extend such methods to solve unpolarized PRD 
problems. 

The GS iterative scheme is obtained by choosing c = old and a = b = new in Equa- 
tion dSJ). This gives 

J*,i = ^!f +neW + / //',-V,'.;; SS^i dx', (17) 

where SS Xti = p x 5S lx>i and J£f +new is the mean PRD intensity calculated using the 'new' 
values of the source function at grid points 1,2, — 1 and the 'old' values at points 
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i, i + 1, • • - , N. The line source function corrections are now given by 

8S ixA - (1 - e) / gl xl h x iAiPx>8S lxl A<lx' 



= {l-e)J^ + ^ + eB-S^. (18) 

To compute 5Si x j we can apply either the FBF or the CRD-CS method discussed in § 3.1. 
The SOR iterative scheme is obtained by doing the corrections as follows : 

ASgf = ljSSS* , (19) 

where u is a parameter with an optimum value between 1 and 2 which can be found using 
the method discussed in § 2.4 of TF95. The Optimum value of u is the one that leads to the 
highest rate of convergence. We find that the SOR method cannot be combined with the 
CRD-CS method for type II redistribution, while it works fine with the FBF method. The 
reason could be the way the wings are handled in the CRD-CS method. 



3.2.1. The standard GS and SOR techniques 



It is worth noting that the GS iterative scheme is twice faster compared to the Jacobi 
scheme. A factor two of additi onal improvement ca n be achieved by implementing the sym- 
metric GS scheme (see TF95; iTrujillo Buenol 120031 ) . To explain the symmetric GS scheme, 
we first recall the GS scheme briefly, which is explained in greater detail in TF95. 

Following TF95, we consider two distinct parts: a incoming and outgoing. 



1. Incoming part (/i < 0) : One starts from the upper boundary (i = 1 with % being 
the depth index), and determines the intensity of the incoming rays (// < 0) at all depths 
using the short-characteristics formal solver. Thus, at the end of the incoming section, one 
has calculated the incoming contribution to the mean PRD intensity, J x a{^ < 0), at all the 
depth points (i — 1, • • • N). 

2. Outgoing part (/x > 0) : One now starts from the lower boundary (i = N). Given that 
at this point the intensity is known, one can easily compute the total mean PRD intensity 
J x ,n- We then use it to calculate SSi X) n and thereby the new source function S^™ at the lower 
boundary. Now for the next depth point i = N — 1, to calculate Jjv-i using Equation ((7|) GS 
uses S^, S^'at.x and 5'° 1 ^ r _ 2 . Then we compute the outgoing contribution to J^jv-i^ > 0). 
However, note that the incoming contribution to J x .n-\{^ < 0) was calculated with the old 
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values of the source function, namely S° t %, S° and S°r^_ 2 . Therefore, to calculate the 
actual J° l< ^i W we add the following correction : 

AJ^x = J g^Af^dx', (20) 

where 

i r° 

A ^iv-i = ^S X , N J ^ x>N (fi < 0)d// . (21) 

Here "in" stands for the incoming pass (see below). Once the actual J^f^ w is found 
one computes SS^n-i and the new source function S l ° e J y_ 1 . Because S^jv-i is available now, 
before going to the next depth point the following correction should be added to the intensity 
I^ >0 (N-1): 

&I%>o(N - 1) = **, N -i{fJ. > 0)5S XtN ^. (22) 

The above procedure is then repeated for subsequent depth points. Clearly, unlike the Jacobi 
iterative scheme, the GS iterative method requires specific ordering of loops in the formal 
solver. The outermost loop is over the directions (first the incoming and then the outgoing 
rays). The next loop is over the spatial points, followed by the loop over different \n\ points. 
The innermost loop is over the frequencies. 



3.2.2. The Symmetric GS and SOR techniques 



The symmetric GS iterative scheme discussed by lTrujillo Buenol (120031 ) is obtained by 
introducing one more outer loop, which first does the GS iteration starting with the incoming 
ray (which we call incoming pass), and then the GS iteration starting with the outgoing ray 
(which we call outgoing pass). In the incoming pass all the GS steps are exactly the same 
as those described above. In the case of the outgoing pass, we again consider two parts to 
describe the symmetric GS case, namely the outgoing and the incoming parts. Let us clarify 
these parts : 



3. Outgoing part of the outgoing pass (// > 0) : We start from the lower boundary 
(i = N) and compute the outgoing contribution to the J x ,i{^ > 0) quantity at all depth 
points using the source function computed newly from the incoming pass. 



4. Incoming part of the outgoing pass (fi < 0) : We now start from the upper 
boundary. At the upper boundary (i — 1), we can compute the new source function S^™, as 
the intensity is known. For the next depth point i = 2, we compute the intensity using S^, 
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S$$, and S° l l Note that the so-called S° l f for i > 2 are the new source functions obtained 
from the incoming pass (see above). Thus, we can now compute the incoming contribution 
to the J x ^{^ < 0). However, the outgoing contribution to J x ^{^ > 0) was calculated with 
the old values of the source function, namely, S° l f, S° l 2, and S^. Therefore to calculate the 
actual J x 2 new , one has to add the following correction 



^Ajr 2 d*' 



where 



AJ° ut 



1 



'x,2 



+ 1 



5S xA / > 0)du. 



(23) 



(24) 



In the above equations "out" denote the outgoing pass. Once the actual J° l $ +new is found 
we can now compute the new source function S^ e 2 w . Since S^™ is available, before going to 
the next depth point the following correction should be added to the intensity I XA1 <o(2) : 



AJJ*o(2) = <M// < 0)5S Xi2 . 
The above procedure is then repeated for subsequent depth points. 



(25) 



This scheme together with the incoming pass (1 and 2) and outgoing pass (3 and 4) is 
nothing but the symmetric GS iterative scheme (hereafter SYM-GS). Thus, each call to the 
formal solver produces as an output two truly GS iterations. Clearly, the incoming pass has 
a convergence rate equivalent to that of a lower triangular approximate operator method 
and the outgoing pass has a convergence rate equivalent to that of an upper triangular 
approximate operator method (see TF95). This s ymmetric GS scheme can be extended to 
SOR also, which is then called Symmetric SOR (jTrujillo Buenol 120031 ). The advantage of 
SSOR is that i t is less sensitive to t he choice of the optimum value of u as compared to SOR 
sec Fig. 1 of iTruiillo Buenol l2003l l Furthermore, unlike S OR, the SSOR m ethod can be 
combined wit h standard acceleration techniques like Ng (s ee lAuerlll987l . Il99ll ) or orthomin's 
acceleration (jVinsomelll976l ; iKlein et al.lll989l ; lAuerlll99lh . 



4. The true error of the numerical solutions 



Following lAuer et al.l (119941 ) we define three quantities that characterize any iterative 
scheme, namely (1) the maximum relative change R c , (2) the maximum relative convergence 
error C e , and (3) the maximum relative true error T e . For a given level of grid resolution g 
at the nth iterative stage these three quantities are defined as follows : 

\Si x (n,g) - S lx (n - l,g)\ 



Rc(n,g) 



max 

T,X 



Six(n,g) 



(26) 
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T e {n,g) 



C e (n,g) 



max 

T,X 



max 

T,X 



\Six(n,g) - Si x (oo,g)\ 
Six(oo,g) 

Six(n,g) - ^(00,00)1 
Si x (oo, 00) 



(27) 



(28) 



In the above equations (n = 00, g) indicates that one is dealing with the fully converged 
solution on a grid resolution level g, while (n = 00, g = 00) indicates the true solution on a 
grid of infinite resolution. T e (oo,g) is nothing but the truncation error corresponding to a 
grid of finite resolution level g, and thus it determines the accuracy of the converged solution 
in that grid. 

In this paper we find the fully converged solution on a given grid resolution level g, by 
iterating until R c < 1CT 10 . Beyond this value R c does not decrease any further, but simply 
fluctuates around it. The true solution required to calculate the true error is found by using 
a grid which is twice finer compared to the grid on which we seek the true error. 

Following TF95, in this paper we use the true error to determine the convergence prop- 
erties of the iterative schemes. Here we show that the true error not only depends on the 
resolution of the spatial grid and the accuracy of the formal solver, but also on the choice of 
the redistribution function. In the following subsections we discuss the true error separately 
for the -Rijiju.aa functions and for cases with a linear combination of i?n ( AA an d -Riii.aa- 



We recall that physically this case represents an atom with two sharp upper and lower 
levels. Thus, the line is infinitely sharp in the rest frame of the atom. In the laboratory frame 
it is broadened by the Doppler effect. This idealized case can hardly be applied to interpret 
any spectral lines, nevertheless it is an interesting academic case to study as it allows to 
examine the effects of pure Doppler redistribution by a Maxwellian velocity distribution. 

Figure [1] shows the convergence properties of different iterative schemes discussed in 
the previous section, applied here to type I redistribution. For all computations presented 
in this paper we consider the case of a semi-infinite atmosphere with the lower boundary 
condition I Xfl (r = T) = B, and upper boundary condition I Xfl (r = 0) = 0. The depth 
grid is constructed using the relation r = exp(— Z), where Z = z/H (with H the scale 
height), and Z the height in units of H. We choose a uniform spacing of AZ. For all 
the figures presented in this paper we have chosen AZ = 0.25 (which corresponds to 9 
points per decade). A Gaussian quadrature with 3 inclinations [0 < /i < 1], and an equally 
spaced frequency grid with 41 points and a spacing of 0.25 Doppler widths are used. Note 



4.1. Pure Doppler redistribution - type I redistribution 
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that a frequency bandwidth of < x < 10 is more than sufficient, as we are considering 
a pure Doppler redistribution (with zero damping). The collisional destruction probability 
e = 10 -4 . The Plank function B is set to unity. Unless stated otherwise we set the continuum 
parameter r to zero. From the left panel of Fig. [TJ we see that the convergence behavior of 
the different iterative schemes are exactly the same as that discussed in TF95. We note that 
for type I redistribution we obtain a true error of 2.7 x 10~ 3 , while for the corresponding 
coherent scattering and CRD (with damping parameter a = 0) cases we get a true error of 
3.5 x 10~ 3 , and 4.3 x 10~ 3 , respectively. 

4.2. Doppler, Natural and Collisional Broadening - type III redistribution 

Physically this case represents a resonance line with its upper level both radiatively and 
collisionally broadened. Collisions are so frequent that there is CRD in the rest frame of the 
atom. 

Figure [2] shows the convergence properties of the different iterative schemes, applied 
here to type III redistribution. Model parameters are the same as those for type I redistri- 
bution, but now the damping parameter a = 1CT 3 . The angular and depth grids used for 
the computation are exactly the same as those used for type I redistribution. However, a 
non-uniform frequency grid that extends up to 1000 Doppler widths from the line center is 
used (as now a ^ 0). We note that for type III redistribution we obtain T e = 2.3 x 10~ 3 , 
which is nearly the same true error as that obtained for the corresponding CRD case (with 
a = 10 -3 ). This is expected, as it is well known that i2ni,AAj i n t ne res t frame of the atom 
behaves like CRD. 



4.3. Doppler and Natural Broadening - type II redistribution 

Physically type II redistribution represents the case of a line with an infinitely sharp 
lower level and an upper level broadened by radiative decay only. In the rest frame of the 
atom the absorption profile is a Lorentzian and the scattering is completely coherent. This 
type of scattering problem is essential to model strong resonance lines, formed in low density 
media. 

Figure E] shows the convergence properties of different iterative schemes, applied here to 
the type II redistribution problem. The model parameters and the various grids used for the 
computation are the same as those used in § 4.2 for type III redistribution. We point out 
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that the SSOR method works well when combined with the FBF technique, while it doesn't 
work properly when combined with the CRD-CS method. The reason is probably due to 
the way the wings are handled in the CRD-CS method. For type II redistribution we obtain 
T e = 0.12, which is pretty a high value compared to that obtained in the type I, III and 
CRD cases. It is well known that one needs a much more refined frequency grid for -Rh,aa 
than for the other redistribution functions because the asymptotic large scale beha vior of 



the t ransfer equation for i?n,AA is like a space and frequency diffusion equation (see iFrisch 



19881 ). However, we checked that use of a frequency grid even finer than the non-uniform 
frequency grid mentioned above does not change the T e value quoted above. Such a high 
value of T e could be due to the fact that -Ri^aa has coherent peaks in the wings, while other 
functions do not have coh erent peaks. Furthermore, in the case of -Rh,aa the wings cannot 



be easily thermalized (see IFrisch! 11980). It is worth to note that very far in the wings only 
diffusion in space remains. Such a regime is encountered only in pure -Ri^aa problems. The 
presence of a background continuum or of some elastic collisions will hide this very far wing 
regime and thereby decreases T e (see below). 

We made a detailed study of the true error for the -Rh,aa redistribution function case 
using the SYM-GS iterative method. Figure HJ shows the true error for different e values. 
Note that the true error decreases when the non-LTE parameter e increases (i.e., when the 
number of scattering events decreases). In Table [1] we present the true error for different 
resolutions of the depth grid. As expected, the true error decreases as the grid resolution 
increases. 

Figure [5] shows the behavior of the true error for type II redistribution when a back- 
ground continuum is included. Clearly, the addition of the continuum decreases the true 
error substantially, as the wings can then be thermalized. Note that even with an opacity 
ratio r as small as 10~ 12 the true error decreases to nearly 3.5 x 10~ 3 , from T e = 0.12 for 
the pure line case. Since in practical problems a background continuum is always present, 
we can conclude that the true error of the numerical methods based on operator splitting 
for type II redistribution can be made significantly small. For example, it is 0.2% when 
r = 10" 4 and AZ = 0.25. 



4.4. Linear combination of -Rhaa and -Rihaa 



We now consider a more realistic case chara cterized by th e following weighted combi- 
nation of type II and type III redistribution (e.g., IStenflolll994j ) : 



Raa(x, x') = 7-Rh,aa(z, x') + (1 - 7)-Riii,aa(^, x'), 



(29) 
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where 7 = 1/(1 + Pe/Tr), with Te the elastic collisional rate and Tr the radiative rate. 

Figure M shows the behavior of the true error for different choices of the elastic collision 
parameter Fe/Tr. The true error corresponding to Te/Fr = is nothing but that corre- 
sponding to the pure -Rh,aa case, which shows the largest value for the truncation error. 
Introducing a small mix of type III redistribution through the contribution of elastic colli- 
sions results in a decrease of the true error. Already for Te/Fr = 0.1, the true error is nearly 
the same as that corresponding to the CRD case. This again shows that the coherent peaks 
of -Ri^aa are responsible for a large truncation error in the case of pure type II redistribution 
without any background continuum. 



5. Polarized PRD radiative transfer equation 

In this paper we restrict ourselves to situations where the radiation field is axially 
symmetric. This condition is satisfied only for one-dimensional plane-parallel or spherical 
atmospheres with either no magnetic field, or a micro-turbulent and isotropic field, or a micro- 
structured magnetic field with a fixed inclination and a random azimuth. Here we consider 
the case of a plane-parallel atmosphere with zero magnetic BeldQ An ax ially symmetr ic 



polar ized radiation field is described by the Stokes parameters I and Q (see I Chandr asekhar 



19501 ). where / denotes the intensity and Q the linear polarization (i.e., the difference between 
the intensity components parallel and perpendicular to a given reference direction in the plane 
perpendicular to the direction of the ray under consideration). In this paper, the positive Q 
direction is defined in the plane containing the direction of the ray and the vertical Z-axis. 
The one- dimensional transfer equation for the Stokes vector components I XfJit j = (I, Q) for 
j = 0,1 is given by 

frtxnj ( r ) = Ixnj (t) - S Xfl ,j (r) . (30) 
The source vector components S XfM j = (5^„, S^J for j = 0, 1 are of the form 

_ + rBUj 



1 We remark that a microturbulent magnetic field can be taken into account by simply replacing 
W2(Ji, J u ) (see § 5.1. for its definition) by H2 Wvt{Ji, Ju), where H2 is the so-called Hanle depolariza- 
tion factor. H2 is unity when the magnetic strength is zero. The explicit form of H2 for a n isot ropic 
magnetic field and a horizontal magne tic field with random azimuth can be found in IStenflol (jl994h and 



Landi Degl'Innocenti fc Landolfil (|2004l ). 
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wher e Uj = (1,0) for j = 0, 1, and the line source vector components Si XfJi j are given by 



e.g., 



Rees fc Saliballl982h 



/+oo ^ r + 1 

dx'~ J d/j,' 

i 

x ^ [R(x, a/; //, [i')] jjf 4yj'. 
i'=o 



(32) 



In the above equation [R( ')] ••, are th e elements of the scattering redistribution ma - 

trix R(x, x'; /i, //) for the non- magnetic case (IRees fc Saliballl982l : iDomke fc Hubenylll988l ). 
In the following subsections we first discuss the redistribution m atrix fo r the n on-magnetic 
case, and then present the decomposition technique proposed by iFrischl (120071 ). This is be- 
cause the iterative algorithms given in this paper are based on the ensuing equations deduced 
in 8 5.2. 



5.1. Redistribution matrix 



A hybrid approximation to R(x, x'; /i, /i') was introduced by IRees fc Salibal (119821 ) : 

R(x, x'- fi, fjf) = (1 - 6)<&,P(m, »'), (33) 



where the pha se matrix P(/x, /x') is given by (e.g., lLandi Degrinnocenti fc Landolfil I2004J ; 
Bommierlll997l ) 

P(/i,^)= W K (Ji,Ju)P%(»,v')- (34) 

K=0,2 

The coefficient Wq(Ji, J u ) = 1, with Ji and J u being the total angular momentum quantum 
numbers of the lower and upper levels, respectively. The coefficient W 2 (Ji, J u ) characterizes 
the maximum linear polarization that can be produced in the line. In the case of a normal 
Zeeman triplet (J; = 0, J u = 1), W2(Ji, J u ) = 1, and P(/i, //) = Pr(/i,/i') is the so-called 
Rayleigh phase matrix. Even though the figures of this paper correspond to the case of a 
normal Zeeman triplet, we present the equations for the more general case of arbitrary values 
of Wk( Ji, Ji,)- The Rayleigh p hase matrix multipolar components P^; (/i, //) are given by 



sec 



Landi Degl'Innocentilll984l . written here for the azimuthally symmetric case) 



[p£(a*,aO]^ = % K (j,e)f K (f,e'), 



(35) 



where j,j' = 0,1. The notation 7q (j, 6) was introduced by IFrischl (120071 . see her Equa- 
tion (28)), where for each K, Q takes values between —K to +K in steps of unity. These 
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quantities are related to th e irreducible tensors for polarimetry Tq{j, f2) i ntroduced b y 
Landi DeglTnnocentil ( 119841 ). where f2 = (9,x) denote the ray direction (see iFrischl 120071 ). 
Since here we are dealing with the azimut hally symmetric case, the relevant qua ntities cor- 
responding to Q = are (see Table 5.6 of lLandi Degl'Innocenti &: Landoml 120041 ) 



7?(O,0) = 1; t 2 (0,#) 
7?(1,0) = O; 7?(1,0) 



1 



2V2 
3 



2^ 



(3^ 2 -l), 
(l-/i 2 ). 



(36) 



In Equation (13"3"|) . e = ^/(I^ + T R ) with 1^ being the inelastic collisional rate and T R 
the radiative rate. However, Equation fl33|) is only an approximate form of the redistribution 
matrix as it does not take into account the effect of elastic (Te) and depolarizing (D^) 
collisional rates. The first quantum mechanical calculation of the redistribution matrix for 
th e resonance po l arizat ion, taking into account the effect of elastic colli si ons was performed 
bv lOmont et al.l ( 119721 ) . Starting from the work of lOmont et al.l (119721 ) . iDomke fc Hubeny 
f)l988h derived a tractable analytic expression of the red istribution matrix. It is worth to 



note that the redistribution matrix that was derived by IDomke fc Hubeny! ( 11988 



is very 



Hummer 



genera l, namely, it depends on the angle-dependent red istribution fun c tions of 
( 1962 ). However for computational simplicity, following iRees fc Salibal (1l982[ ). iNagendra 
( 119941 . see also Faurobert-Scholl 1992) used the angle-av eraged version of the Domke-Hubeny 
(DH) redistribution matrix. Following iBommierl ( 119971 ) we write this redistribution matrix 
as follows : 



where the branching ratios a and are given by 

r R 



3 



a 



r R + ri + r E : 
r R 



(37) 



Note that 



r R + r : + DW 

0, and also that the factor (1 — e) is contained in the branching ratios. 



(39) 



It is worth to clarify certain important points related to these branching ratios (see also 
Nagendralll994l ). In astrophysics one expects that the branching ratios add up to unity. How- 
ever, from Equation ( 137|) we see that the branching ratios add up to give [a + 
fr K \ which for K = is nothing but (1 - e) and for K = 2 is (1 - e)/ [l + 5 (2) (1 - e)] with 
S (2) = £)(2)/r R . We note that these are indeed the factors that appear in the line source 
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function expressions for Stokes I and Stokes Q (namely S% and S% or p® and p%) in the CRD 



case formulated (see iTrujillo Bueno &: Manso Sainzl Il999t lLandi Degl'Innocenti fc Land olfi 



2004). It is important to note that some authors (e.g. 
19941 ). write the factor (1 — e) before the second term of Equation 
branching ratios a and ^ K > by (1 — e), namely 



Faurob ert- S choll 



1992 



Nagendra 



and renormalize the 



a 



i old 



old 



" l-e ~ r R + r! + r E ' 
P {K) _ r R + r t 

l-e T K + T 1 + DW 



(40) 



The approximate form of R(x, x'\ p, p') gi v en in Equa tion (1331) was used in pl ane-parallel 
polarized radiative transfer by lRees Salibal (1l982h and iFaurobertl (119871 . Il988l ). with g x \,. 
These authors used Feautrier's (19 64) method to solve the pola rized transfer equation. A 
discrete space method wa s used by Nagendral (119861 . Il988l . Il989l ) for the same problem but 



in spherical atmospheres. iMcKennal (119841 ) used an integral equation approach to solve the 
same problem with a linear combination of g xx i and g xx i (see Equation ([29])). The DH redis- 
tribution matrix given in Equation (1571) was use d in plane-parallel polarized radiati ve transfer 
computation s by iFaurobert-SchollI ( 119921 . Il993l ) . The same problem was solved by iNagendra 
( 119941 . I19951 ) but in spherical atmospheres. As already mentioned in the introduction such 
methods are computationally expensive. 

A Jacobi based ALI method to solve the polarized radiative transfer e quation with the 
hybrid approximation for the redistribution matrix and g l xx , was developed by 
( 



Paletou &; Faurobert-Scholl 



1997). They extended the CRD-CS method of PA95 to scattering polarization. ITrujillo Bueno &: Manso S 



(119991 ) generalized the symmetric GS and SOR methods of TF95 to CRD polarized radiative 



transfer, with the relevant equations formulated within the framework of the quantum theory 
of spe ctral line polarization described in the monograph by lLandi DeglTnnocenti fc Landolfi 
(120041 ) . In this paper we generalize these symmetric GS and SOR methods to solve the 
above-mentioned PRD problem. We consider both, the hybrid approximation and the DH 
redistribution matrix. In the next section we present the Jacobi, GS and SOR iterative 
methods to solve polarized radiative transfer problems with the DH redistribution matrix. 



5.2. Decomposition in the irreducible basis 

From Equations (|3T!) and (132]) we see that unlike the unpolarized case, the line source 
vector components now depend not only on the frequency x but also on the orientation /x 
of the radiation beams. In the case of CRD the line source vector components depend only 
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on At and are independent of x. To reduce the computational cost, iFaurobert-Scholl et al. 
(119971. see also P aletou & Faurobert-Scholl 1997) used a factorized form of P(/i, //) given 
by llvanovl ( 119951 ). which allowed them to transform or reduce the polarized CRD transfer 
equation to a 2 x 2 basis wherein the source vector components are independent of \x. To this 
2x2 matrix transfer equation they applied a Jacobi iterative scheme to solve the problem. 

The factorization of the Rayleigh phase matrix int o a product of two 2x2 matrices 
that depend separately on \i and // is not unique (see iFrischl 120071 ). Such a factorization 
come s out naturally if one uses the Tn(j, fi) irreducible tensors (see lLandi DeglTnnocenti 
19841 ) to d erive the Ray leigh phase matrix (see Equation (1351) ). Using the irreducible tensors 
Tn{ji H). IFrischl (120071 ) provided a simple way of transforming or reducing the Stokes vector 
components to irreducible tensors in the case of the Hanle effect regime. Such a transforma- 
tion is referred to as the "decomposition" of the Stokes vector components. We note that 
such a deco mposition comes out naturally in the de nsity matrix theory of spectral line polar- 
ization (see lLandi DeglTnnocenti fc Landolfill2004l ). Furthermore, it is well known that the 



density matrix and the scattering formalisms (that we adopt in this paper) are equivalent for 
a two-level atom without lower-level polarization and stimulated em ission. However, in this 
paper we use the decomposition technique proposed by IFrischl (120071 ). but applied here to the 
axially symmetric case. For clarity, we p resent some im portant steps of this decomposition. 
For more details the reader is referred to IFrischl (120071 ). 



For the azi muthally symmetric case the Stokes vector component decomposition given 
by IFrischl ( 120071 ) takes the following form (in the notations used in this paper) : 



K=0,2 



(41) 



K 



with similar equations relating Uj, S Xfl j and Si Xfl j to U , 



(S x )q and (Si x )q, respectively. 



Note that U$ 



0. The quantities (I xl i)o and (S x )q are called the irreducible 
tensor components of the Stokes and the source vector components, respectively. 



1 and C/ 2 



Substituting Equations ( 155]) . (J57|) and (l4"Tj) in Equations (1501 - (152]) . it can be shown that 
(I X fi)o satisfies a transfer equation similar to Equation ( 15U|) but with I Xfl j and S XfM j replaced 
by (Ix/j,)o and (S xij )q , respectively. Furthermore, (S x ^)q is given by Equation (|3T|) but with 
Si XfJi j and Uj replaced by (S^jf and Uq , respectively. The irreducible components of the 
line source vector are now given by 



eBUK + W K (J h J u )(J x )X, 



(42) 



where 



(Jx)o 



dx' 
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>x>)q- 



In the above equation the angle integrated irreducible tensor is given by 



(Jx')o 



K 



K'=Q,2 



where 



(43) 



(44) 



(45) 



i'=o 

Clearly, the advantage of this decomposition is that the irreducible tensor components of the 
source vectors are now independent of the orientation \i of the radiation beam. 



Following iFrischl (120071 ) we now introduce the 2-component Stokes and source vectors 
in the irreducible basis 



21 T 



[(^)o' (^)o] 



21 T 



(46) 



In the above vector notation, the transfer equation in the irreducible basis can be written as 



1-xiSj) - S x (t), 



(47) 



where S x is given by Equation (13~T1) . but with Si xli) j and Uj replaced by Si x and U, respec- 
tively Here U = (1,0) T , and 

S lx = eBU + WJ x . (48) 



In the above equation 



where 



J, 



^xx'Jx' dx' 



N^^g^aE + g™, (B - aE) 
and the 2-component mean intensity vector 



1 f +l 

- J ^ *(//) Z^d//. 



(49) 
(50) 

(51) 



In Equations (l5"Uj) and (ISTi) . E is a 2 x 2 identity matrix, while W = diag[Wo, and 
B = diag[/3 (0) ,/3 (2) ] are 2 x 2 matrices. Note that since the matrix B is diagonal, the matrix 
N xx ' is also diagonal. The expli cit form of th e 2x2 matrix \I/ formed by the elements ^q K ' 
can be found in Appendix A of IFrischl (120071 ). In the next section we apply the Jacobi, GS 
and SOR iterative schemes to Equations ()47l)-(j5T!). 
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6. Iterative methods for polarized PRD radiative transfer 

The formal solution of Equation ( |4"T|) is given by Equation (JSJ), but with 1^, A^, 5* x 
and T x ^ replaced by X XA1 , A Xfl , S x , and T x/J , respectively. Here T XfJi is the transmitted 
2-component Stokes vector due to the incident radiation at the boundaries, and is a 
2N x 2N operator. For given depth indices A XfJi} a> is a 2 x 2 block. We use again the 
short-characteristics method as the formal solver, but now applied to the vector transfer 
equation fl3J~ 



As in Equation ((§]), we now write the 2-component mean intensity vector as 

3" x,i = A.x,il&x,l "I" ' ' ' "I" A-x,ii—l&x,i—l 

+A Xjii S Xii + A Xjii+ iS c Xii+1 + ■•• 

+A Xj jiv«S^ 7V + T X) i , (52) 

where T x i is given by Equation ( 15 ip but with I x ^ replaced by T XfJj . The 2x2 matrix A Xtii i 
is given by 



1 f +1 

A x ,n' = - J A^ji' d/i. 



(53) 



Note that unlike the unpolarized case (see Equation (jSJ)), at each depth we now have to 
perform a matrix operation involving a 2 x 2 matrix and a 2 column vector. In the following 
subsections we successively present the Jacobi, GS and SOR iterative schemes. 



6.1. Jacobi Iterative scheme 

It is straightforward to generalize the Jacobi scheme discussed in § 3.1 to the scattering 
polarization case. With this iterative scheme, the equations for the 2-component line source 
vector corrections are given by 

/+oo 
NrWvA.,.'.,, 5S lx , yi dx' = TZ X)i , (54) 
-oo 

where 

n x , = eBU + WJ^-S? x %. (55) 

As discussed in § 3.1. the system of linear equation f}54"]) can be solved by the FBF method, 
namely 

A5S = n, (56) 
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where at each depth point i, TZ is a vector of length 2iV x , and the matrix A. is of dimension 
2N X x 2A^.. For a given depth point i, and given frequencies x,x', A. is a 2 x 2 block denoted 
by ^4 2 , and given by the expression 



3mn E N mn J) n A r 



m, n 



AL. 



(57) 



Clearly, the above FBF method is numerically much more expensive compared to the unpo- 
larized the size of the matrix A . is now twice larger. The FBF method was actually 

generalized by ISampoorna et al.l ( 120081 ) for t he weak field regim e of the Hanle effect (with 
PRD). This method has also been used by iFrisch et al.l (120091 ) for the Hanle effect of a 
turbulent magnetic fields with a finite correlation length (and with CRD). 

As already discussed in § 3.1. for the unpolarized case, the CRD-CS method is compu- 
tationally less expensive, but is as robust as the FBF method. This CRD-CS meth o d was 
extended to scattering polarization (non-magnetic) by iPaletou &: Faurobert-Scholll (119971 ) 
for the hybrid approxima t ion th at uses g l xx ,. This method was later extended to the Hanle 
effect cas e by iFluri et al.l (120031 ) . They used the weak field Hanle redistribution matrix of 
Bommierl ( 119971 ). the so-called approximation-Ill, which reduces to the DH redistribution 
matrix for the zero magnetic field ca se. In the following subsection we briefly recall the 
CRD-CS method of IFluri et al.l (120031 ). applied here to non-magnetic case. The main differ- 
en ce with the e quations presented in this paper is that we use the decomposition technique 
of iFrischl ( 120071 ). whil e IFluri et al.l (12003!) us e the traditional Fourier-azimuthal expansion 
technique discussed in iNagendra et al.l ( 119981 ) . 



6.1.1. CRD-CS or core-wing method for the DH redistribution matrix 

As discussed in § 3.1.1 and 3.1.2, in Equation (151)) we approximate g\\i by Equation (fl4l) 
and g™, by <p x i in the line core (x < x c ), but set it to zero in the wings (x > x c ). This gives 
(using Equation (|50|) ) 

(58) 



5S 



lx,i 



1 - a x a p x W A X)i 



Following IFluri et al.l ( 120031 . see also Paletou & Faurobert-Scholl 1997; Nagendra et al. 1999), 
one can show that 



AT, 



where 



E- 



+x c 



dx(j) x p x A xU ■ WB 



x Px Ajj^j 7Zx,i dx. 



(59) 



(60) 
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Note from Equation fl58|) that in the core (when a x = 0) the denominator reduces to unity 
(i.e., we are left with a simple summation), while in the wings the term multiplying (1 — a x ) 
appears as a frequency independent quantity. 



6.2. GS and SOR iterative schemes 



These radiative transfer methods were developed bv lTruiillo Bueno &: Manso Sainzl (11 999) 



for the non-magnetic and micro-turbulent field CRD cases and by lManso Sainz &: Trujillo Bueno 



(119991 ) for the CRD case of a deterministic magnetic field in the Hanle effect regime. Here 
we present the generalization of these methods to the PRD problem of resonance line polar- 
ization in the absence or in the presence of a weak magnetic field which does not break the 
axial symmetry of the problem (e.g., a microturbulent field). 

The GS iterative scheme is obtained by choosing c = old and a = b = new in Equa- 
tion (IB^|) . which gives 

U x,i — *T° 1 i i +new + A.x,u 5S Xt i , (61) 

where J° x ^ +new is the 2-component mean intensity vector calculated using 'new' values of 
the source vector S Xji at grid points 1, 2, • • • , % — 1, and the 'old' values at i, i + 1, • • • , N. 
The line source vector corrections are given by Equation (IM|) . but with 

n x , l = eBU + WJ^-S^. (62) 

The line source vector corrections can be computed using either the FBF or the CRD-CS 
methods discussed in § 6.1. Note that the GS as well as the SYM-GS iterative algorithms 
discussed in § 3.2 can be extended straightforwardly to the polarized case. The only difference 
is that we are now dealing with a 2-component source vector, with intensity as well as mean 
intensity vectors, and with an approximate lambda operator which is now a 2 x 2 block for 
any given depth and frequency. For example, in the several correction terms that one needs 
to consider in a SYM-GS algorithm, namely Equations (|20l)- (p5j) . we have to replace simply 
the unpolarized quantities J Xj i, J x ^ and S x> i by the polarized 2-component vectors J~ x ^, 3 
and «S X] j, respectively, to be able to apply those equations to the polarized case. Therefore 
unlike in the unpolarized case we now have to do several matrix manipulations (see, e.g. 
Equations ( 156]) and (l58|) ). 



X.l 



The SOR iterative scheme is obtained by doing the corrections as follows : 

6S^ = u6SZ- (63) 



As already noted, all the three iterative schemes (Jacobi, GS and SOR) involve matrix 
operations for the computation of the line source vector corrections (see Equations (|56|) 
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and fl58l) ). A smart strategy to a void such matrix computations , and thereby speed up 
the iterative methods, was given by lTrujillo Bueno fc Manso Sainzl ( 119991 ). To describe this 
strategy in some detail, we now write the 2x2 approximate operator A x u as follows : 



A 00 . A 

T' 11 



*-X,ll 



,00 

x,n 

.20 



02 

x,ii 

, 22 



A 20 - A 2 

XM XM 



(64) 



where A x ^' are given by (see Equation (1531 ) 



A 



x,ii 



+ 1 



(65) 



As shown by lTrujillo Bueno fc Manso Sainzl (119991 ) for the CRD problem (see their Fig. 3), 
we also find for our PRD problem that |A°^| > |A 22 J > |A° 2 J = |A^|. We illustrate 
this fact in Fig. [7J where the elements of the monochromatic lambda matrix for a semi- 
infinite model atmosphere is plotted versus the line center optical depth for two different 
frequencies x = (panel a) and x = 5 (panel b). Clearly, for x = the elements of A x ii 
ar e nearly identical to that of the corres ponding CRD case (compare our Fig. [TJi with Fig. 3 
of iTrujillo Bueno &: Manso Sainzl Il999l ) . As the frequency increases toward the line wing 
the entire curve corresponding to all the elements of this matrix shifts toward higher optical 
depths (see Fig. [7b). In other words, the depth at which A^' reaches unity, 0.7 and zero 
for (K,K') = (0,0), (2,2) and (0,2), respectively, shifts toward larger optical depth. This 
behaviour can be understood by looking at the explicit form of A x ^. In the case of a short- 
characteristic formal solver, A x/ f^ ^^(a*), where O is the depth point of interest (here 
O = i). The quantity ty Xt o(fJ>) is given by 

(Arp - At m )wi + w 2 



At p At m 



(66) 



with wo = 1 — exp(— Atm), w\ = w — Atm exp(— Atm) an d w 2 = 2wi — Arfj exp(— Atm)- 
Clearly, when both Atm and Arp tend to infinity, ^ X) o{^) — > 1, and then it is easy to see 
from Equations (|65|) . fj4"5]) and f l36]) that A x ^' saturate to their respective values mentioned 
above. As x increases the optical depth at which ty X} o(fi) — > 1 shifts to larger optical depths, 
and hence the observed behaviour. Thus, saturation is reached when ^f x ,o(l^) ~^ 1- This is 
perhaps equivalent t o saying that saturation is re ached when the exponential in the kernel 
(see Equation (19) of lFaurobert-Scholl et al.lll997| ) can be replaced reasonably well by a delta 
function in the grid interval around the optical depth Tj. It is worth noting that as A x ^' 
depends on the optical depth grid, the finer the grid the slightly larger is the depth at which 
saturation is reached. 



Thus, following ITrujillo Bueno fc Manso Sainzl (119991 ) we can also set A 



22 

xM 



A 



02 



A 20 . 

XM 



0, and still obtain a radiative transfer method with a convergence rate that is as 
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go od as that achieved by keeping the f ull A Xt u given in Equations (164"1) and fl65|) (see Fig. 4 
of iTruiillo Bueno fc Manso Sainzlll999h . The use of such a strategy leads to a decoupling of 



the equations for (5Si Xti )Q and (SSi Xji )l as follows (see Equations (15"4"|) . (j5"01) . and (JMjl ) : 

= ^ + ( J* d )g - (5g)g , (67) 

(6S lx>i f = W 2 (J x f)l - (5gJ)g , (68) 

where (iV a!a ./)oo is the first diagonal element of N xx / given in Equation ( 1501) . 

In summary, the (5Sta,i)o correction is computed using a Jacobi, GS or SOR iteration, 
while the (SSi X} i)q correction is formally similar to the classical A-iteration. However, the 
important difference with respect to the classical A-iteration method is that the (I Xfl ^)Q which 
enters the computation of (J x ,i)o ( see Equation flBTl) ) is calculated with the improved (S^q 
value that was obtained in the previous iterative step. 

Finally, we remark the following two important points : (1) As in the case of unpolarized 
transfer, in the polarized case we again find that as long as Fe/Tr < 10 we can approximate 
9 l xl' by 4>x' in the core and neglect it in the wings. But as soon as r E /T R > 10, to get the 
converged solution we have to approximate g^, by <f> x ' throughout the line profile for the 
line source vector correction computation. (2) All the above-mentioned iterative schemes 
are given for the DH redistribution matrix. It is not difficult to apply them for the hybrid 
approximation (see Equation (1331 ). In this case we simply replace N xx > by (1 — e)g xx ,'E and 
in Equations ( 1581) and ( 1591) set a = (1 — e) and B = (1 — e)E. Furthermore, when applying 
CRD-CS for k — I, III we do the same approximations as we did for unpolarized case (see 
§ 3.1.2.). 



7. The true error for (S x )l 

In § 4 we presented a detailed study of the true error for the unpolarized PRD source 
function. In this section we present similar studies for (S x )q. We note that the behaviour 
of T e , C e and R c presented in § 4 for the unpolarized source function is similar for (S x )q. 
Therefore, we consider only (S x )l here. The results are presented in Figs. [8]and[9l 

The definition of T e , C e and R c for (S x )l is also given by Equations (f2"6"|) - (T2"8l . but with 
S[ x replaced by (S x )l. However, (S x )l is a sign changing quantity. Thus, in the denominator 
of Equations (T26T) f[28|) . we need to replace Si x by ((S^ol- Furthermore, it is well known that 
in general (S x )q is two orders of magnitude smaller than (5* x )g. Therefore, we find that in 
our PRD case the maximum relative change R c for (S x )q reaches approximately a minimum 
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value of 10~ 8 and then starts to fluctuate around it. Thus, to find a fully converged solution 
on a given grid resolution level g, we iterate until R c < 10 -8 . We remark that we adopt the 
same method as described in § 4 to find the true error as well as the convergence error for 

As in § 4, here we consider the hybrid approximation with i?i,n,iii,AA redistribution 
functions and the DH redistribution matrix. Again, we consider a semi-infinite atmosphere 
with the lower boundary condition (I x ^)o = B5ko, an d the upper boundary condition 
{Ix[i)o = for K = 0, 2. A depth grid of 9 points per decade (i.e., AZ = 0.25) and a 
Gaussian quadrature with 5 /z-values [0 < fi < 1] are used. The frequency grid used is 
exactly the same as that chosen for the unpolarized case. Other parameters are e = 10 -4 , 
r = 0, B — 1, and a = 10~ 3 for type II and type III redistribution, unless stated otherwise. 
We initialize all the three iterative schemes discussed in this paper by the LTE solution : 
(Si x )q = B 5ko- It is worthwhile to note that if the initial or starting solution is other than 
the above-mentioned LTE solution, then the true error that one obtains for a given grid 
resolution (after reaching the plateau region of the T e curve) remains the same. However, 
the path followed to reach that T e is different. 

We find that the true error for (S x )q is in general larger by one order of magnitude 
compared to that for (S x )q. Fig. [8] shows the T e , C e and R c for (S x )q with the hybrid 
approximation and i2i,n,m,AA redistribution functions. As in the unpolarized case, the true 
error for (S x )q is the largest for the Ru,aa redistribution function case. It is about 23%. 
However, as discussed for the unpolarized case, addition of the continuum or inclusion of 
elastic collisions through the use of the DH redistribution matrix improves the true error 
significantly. For example, the true error for (S x )q is approximately 2.3% for Fe/Tr = 1 
(see the bottom solid line in Fig. [9]) without continuum, and it is 1 % for r E /T R = and 
r = 10~ 4 (figure not shown). 



8. Conclusions 



In this paper we have shown how to solve efficiently and accurately the non-LTE reso- 
nance line formation problem in stellar atmospheres, taking into account PRD effects with 

and without scattering polarization. To this end, we have generalized the Gauss-Seidel (GS) 

and Su cces sive Over- Relaxation (SOR) rad i ative t ransfer methods that lTrujillo Bueno fc Fabiani Bendicho 
( 119951 ) and iTrujillo Bueno fc Manso Sainzl (119991 ) developed for solving unpolarized and po- 
larized CRD problems, respectively. These iterative methods are based on the concept of 
operator splitting. As in the CRD case, we find that these methods are superior to the Jacobi- 
based ALI method. Quantitatively, the symmetric GS method (SYM-GS) is 4 times faster 
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than Jacobi, while the symmetric SOR method (SSOR) is about 10 times faster without the 
need of refining the choice of the cu-parameter. We emphasize that our implementation of 
these highly convergent radiative transfer methods do not require neither the construction 
nor the inversion of any non-local A-operator, so that the computing time per iteration is 
similar to that of the Jacobi method. Therefore, these GS-based methods are suitable also 
for the solution of non-LTE problems in three-dimensional model atmospheres. 

For the unpolarized PRD problem, we have considered the case of pure Doppler re- 
distribution (type I), Doppler, natural and collisionally broadened type III redistribution, 
Doppler and naturally broadened type II redistribution, and a combined case of type II and 
type III redistribution. For the PRD problem of resonance line polarization we have con- 
sidered both the hybrid approximation with angle- averaged type I, II a nd H I redistribution 



functions and the general redistribution matrix of iDomke fc Hubenyl (119881 ) that properly 
takes into account the elastic and depolarizing collisions. The methods we have developed 
here can be used also for solving the resonance line polarization problem in the presence of 
a magnetic field that does not break the axial symmetry of the problem. For the case of a 
weak magnetic field with a given strength, inclination and azimuth at each spatial grid point 



the c orresponding redistribution matrices are substantially more complicated (e.g.. lBommier 



19971 ). but the generalization of the same GS-based methods is straightforward in spite of 



the fact that the number of unknowns is three times larger. 

Finally, we emphasize that the PRD radiative transfer problem we have considered here 
is that of a two-level model atom without the possibility of lower-level polarization, which 
implies that it is assumed that only the emission term of the transfer equation contributes 
to scattering polarization. Fortunately, there are several diagnostically important resonance 
lines for which this two- level atom approximation is probably suitable (e.g., the k line of 
Mg II). In forthcoming papers we will show how the application of the computer programs 
described here allow us to gain physical insight and to make predictions on the Q/I shapes 
produced by PRD effects. 
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Table 1. The true error in the case of type II redistribution for different resolutions of the 
depth grid. For all cases, e = 1CT 4 , a = ICT 3 and r = 0. 



AZ 


Number of points 


T 

± e 




per decade 




0.5 


4.5 


0.2571 


0.25 


9 


0.1244 


0.125 


18 


0.0638 


0.1 


23 


0.0533 


0.04 


28.75 


0.0238 


0.02 


57.5 


0.0130 


0.01 


115 


0.0069 




Fig. 1.— Pure Doppler (type I) redistribution. Convergence properties of the various 
iterative schemes in a semi-infinite isothermal model atmosphere with e = 10~ 4 and a spatial 
grid with 9 points per decade (AZ = 0.25). Left panel : solid line (Jacobi), dotted line (GS), 
dashed line (SYM-GS), dot-dashed line (SSOR, u) opt = 1.6). Right panel: solid line (T e ), 
dotted line (R c ), and dashed line (C e ), are computed using the Jacobi iterative scheme. 
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Fig. 2. — Doppler, natural and collisional (type III) redistribution. Convergence properties of 
the various iterative schemes in a semi-infinite isothermal model atmosphere with e = 10~ 4 , 
a = 10~ 3 and a spatial grid with 9 points per decade (AZ = 0.25). The different line types 
are the same as in Fig. [TJ 




Fig. 3. — Doppler and natural broadening (type II redistribution). Convergence properties of 
the various iterative schemes in a semi-infinite isothermal model atmosphere with e = 10 -4 , 
a = 10~ 3 and a spatial grid with 9 points per decade (AZ = 0.25). The line types and the 
model parameters are exactly the same as in Fig. |2j 
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Fig. 4. — Study of the true error for the type II redistribution function problem. Different 
line types: solid line (e = 10~ 2 ), and dotted line (e = 10~ 4 ). A semi-infinite isothermal 
model atmosphere with a = 1CT 3 and a spatial grid with 9 points per decade (AZ = 0.25) 
are used. 
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Fig. 5. — Study of the true error for the type II redistribution function case. Effect of 
the continuum parameter r. The non-LTE parameter e = 10~ 4 and the depth grid spacing 
AZ = 0.25. Solid line: r = 10" 4 , dotted line: r = 10~ 6 , dashed line: r = 10" 8 , dot-dashed 
line: r = 10" 10 , dash-triple-dotted line: r = 1CT 12 , and long-dashed line: r = 0. Note that 
the dotted and dashed lines merge to give a dot-dashed line. 
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Fig. 6. — Study of the true error for the combined case of type II and III redistribution 
function. Effect of the elastic collision parameter I^e/Fr. The non-LTE parameter e = 1CF 4 , 
the continuum parameter r = and the depth grid spacing AZ = 0.25. The topmost solid 
line: I^e/Tr = 0, dotted line: Te/Fr = 10~ 4 , dashed line: I^e/I^r = 10 -3 , dot-dashed line: 
Te/Fr = 10~ 2 , dash-triple-dotted line: r E /r R = 0.1, long-dashed line: r E /r R = 0.25, and 
the bottom most solid line: r E /r R = 1. 
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Fig. 7. — Variation of the diagonal elements of the monochromatic lambda operator with 
the line center optical depth. A semi-infinite model atmosphere with no continuum (r = 0) 
and damping parameter a = 10~ 3 are used. The solid line corresponds to A° ^, the dotted 
line to A°^j and the dashed line to A^. 
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Fig. 8. — T e , C e and R c for (S x )q. Convergence properties of the various iterative schemes. 
Left panels: solid line (Jacobi), dotted line (GS), dashed line (SYM-GS), dot-dashed line 
(SSOR, w op t — 1.6). Right panels: solid line (T e ), dotted line (R c ), and dashed line (C e ), 
computed using the SYM-GS iterative scheme. For type II and type III redistribution the 
damping parameter a = 1CT 3 . A spatial grid with 9 points per decade (AZ = 0.25) is used. 
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Fig. 9. — True error for the DH redistribution function problem. Effect of the elastic collision 
parameter I^e/Fr. The non-LTE parameter e = 10 -4 , the damping parameter a = 10~ 3 , 
r = (pure line case), = 0.5 Tg, and depth grid spacing AZ = 0.25. The topmost solid 
line: Te/Fr = 0, dotted line: Te/Fr = 10~ 4 , dashed line: Te/Fr = 10~ 3 , dot-dashed line: 
Te/Fr = 10 -2 , dash-triple-dotted line: Te/Fr = 0.1, long-dashed line: Te/Fr = 0.25, and 
the bottom solid line : r E /r R = I. 



